Componentes Principales
Retención de variabilidad
library(FactoMineR)
library(factoextra)
library(DT)
acp <- PCA(
X = datos %>% select(Active:Distressed) %>%
select(-c(Apathetic, Distressed, Frustrated, Bored)),
scale.unit = TRUE,
graph = FALSE
)
as.data.frame(acp$eig) %>%
mutate_if(is.numeric, round, digits = 2) %>%
rownames_to_column(var = "component") %>%
datatable(extensions = 'Buttons', options = list(
dom = 'Bfrtip',
buttons = c('csv', 'excel', 'pdf')),
rownames = FALSE)
fviz_eig(acp, addlabels = TRUE)

Contribución de variables
as.data.frame(acp$var$contrib) %>%
rownames_to_column(var = "variable") %>%
mutate_if(is.numeric, round, digits = 2) %>%
datatable(extensions = 'Buttons', options = list(
dom = 'Bfrtip',
buttons = c('csv', 'excel', 'pdf')),
rownames = FALSE)
library(ggsci)
library(tidytext)
as.data.frame(acp$var$contrib) %>%
rownames_to_column(var = "adjetivo") %>%
pivot_longer(cols = -adjetivo, names_to = "componente",
values_to = "valor") %>%
mutate(componente = as_factor(componente),
adjetivo = reorder_within(adjetivo, valor, componente)) %>%
ggplot(aes(x = adjetivo, y = valor, fill = componente,
color = componente)) +
facet_wrap(~componente, scales = "free") +
geom_col(position = "dodge", alpha = 0.7) +
coord_flip() +
scale_x_reordered() +
scale_color_npg() +
scale_fill_npg() +
labs(color = "", fill = "") +
theme(legend.position = "top")

Correlaciones
as.data.frame(acp$var$cor) %>%
rownames_to_column(var = "variable") %>%
mutate_if(is.numeric, round, digits = 2) %>%
datatable(extensions = 'Buttons', options = list(
dom = 'Bfrtip',
buttons = c('csv', 'excel', 'pdf')),
rownames = FALSE)
as.data.frame(acp$var$cor) %>%
rownames_to_column(var = "variable") %>%
mutate(across(where(is.numeric), round, digits = 2)) %>%
pivot_longer(cols = -variable) %>%
ggplot(aes(x = variable, y = name, fill = value)) +
geom_tile(color = "black") +
geom_text(aes(label = value), color = "black") +
scale_fill_gsea() +
labs(x = "", y = "", fill = "Correlación:") +
theme(legend.position = "top",
axis.text.x = element_text(angle = 45, hjust = 1))

Gráficos de componentes
CP1 vs PC2
fviz_pca_biplot(
acp,
col.ind = datos$Moments,
palette = "npg",
addEllipses = TRUE,
label = "var",
col.var = "black",
repel = TRUE,
legend.title = "Moments:",
ellipse.type = "norm",
ellipse.level = 0.5,
select.var = list(contrib = 10),
geom.ind = "point",
pointshape = 21,
pointsize = 2.5,
fill.ind = datos$Moments
) +
theme_minimal() +
theme(legend.position = "top") +
labs(x = "Component 1", y = "Component 2",
title = "")

CP1 vs PC3
fviz_pca_biplot(
acp,
axes = c(1, 3),
col.ind = datos$Moments,
palette = "npg",
addEllipses = TRUE,
label = "var",
col.var = "black",
repel = TRUE,
legend.title = "Moments:",
ellipse.type = "norm",
ellipse.level = 0.5,
select.var = list(contrib = 10),
geom.ind = "point",
pointshape = 21,
pointsize = 2.5,
fill.ind = datos$Moments
) +
theme_minimal() +
theme(legend.position = "top") +
labs(x = "Component 1", y = "Component 3",
title = "")

CP2 vs PC3
fviz_pca_biplot(
acp,
axes = c(2, 3),
col.ind = datos$Moments,
palette = "npg",
addEllipses = TRUE,
label = "var",
col.var = "black",
repel = TRUE,
legend.title = "Moments:",
ellipse.type = "norm",
ellipse.level = 0.5,
select.var = list(contrib = 10),
geom.ind = "point",
pointshape = 21,
pointsize = 2.5,
fill.ind = datos$Moments
) +
theme_minimal() +
theme(legend.position = "top") +
labs(x = "Component 2", y = "Component 3",
title = "")

CP1 vs CP2 vs CP3 - Estático
cp1 <- acp$ind$coord[, 1]
cp2 <- acp$ind$coord[, 2]
cp3 <- acp$ind$coord[, 3]
colors <- pal_npg("nrc")(4)
colors <- colors[as.numeric(datos$Moments)]
library(scatterplot3d)
s3d <- scatterplot3d(
x = cp1,
y = cp2,
z = cp3,
pch = 19,
color = colors,
grid = TRUE,
xlab = "CP1",
ylab = "CP2",
zlab = "CP3",
box = FALSE
)
legend(s3d$xyz.convert(8, 4, 5.5),
legend = levels(datos$Moments),
col = pal_npg("nrc")(4),
pch = 19, box.col = "white")

CP1 vs CP2 vs CP3 - Interactivo
datos$cp1 <- cp1
datos$cp2 <- cp2
datos$cp3 <- cp3
library(plotly)
fig <- plot_ly(data = datos,
x = ~cp1, y = ~cp2, z = ~cp3,
color = ~Moments,
colors = pal_npg("nrc")(4))
fig <- fig %>% add_markers()
fig
Inferencia
Distribuciones
datos %>%
select(Moments, cp1:cp3) %>%
pivot_longer(cols = -Moments) %>%
ggplot(aes(x = value)) +
facet_wrap(~name) +
geom_density()

Boxplots
datos %>%
select(Moments, cp1:cp3) %>%
pivot_longer(cols = -Moments) %>%
ggplot(aes(x = Moments, y = value, color = name, fill = name)) +
geom_boxplot(alpha = 0.5) +
scale_fill_npg() +
scale_color_npg() +
labs(color = "", fill = "")

Anova CP1
- Resultado análisis de varianza:
anova_cp1 <- aov(cp1 ~ Moments, data = datos)
tidy(anova_cp1)
library(emmeans)
emmeans(anova_cp1, specs = "Moments")
## Moments emmean SE df lower.CL upper.CL
## Habituation -0.989 0.486 281 -1.945 -0.0333
## Treatment 1.164 0.209 281 0.753 1.5757
## Break -0.846 0.243 281 -1.324 -0.3684
## Final -1.913 0.486 281 -2.869 -0.9574
##
## Confidence level used: 0.95
- Diferencias entre medias - Tabla:
tidy(TukeyHSD(anova_cp1))
- Diferencias entre medias - Gráfico:
tidy(TukeyHSD(anova_cp1)) %>%
ggplot(aes(x = contrast, y = estimate, ymin = conf.low, ymax = conf.high)) +
geom_point() +
geom_errorbar(width = 0.2) +
geom_hline(yintercept = 0, linetype = 2, color = "red") +
coord_flip()

par(mfrow = c(2, 2))
plot(anova_cp1)

Anova CP2
- Resultado análisis de varianza:
anova_cp2 <- aov(cp2 ~ Moments, data = datos)
tidy(anova_cp2)
emmeans(anova_cp2, specs = "Moments")
## Moments emmean SE df lower.CL upper.CL
## Habituation -1.009 0.331 281 -1.660 -0.358
## Treatment 0.908 0.142 281 0.628 1.188
## Break -0.657 0.165 281 -0.983 -0.332
## Final -1.264 0.331 281 -1.915 -0.613
##
## Confidence level used: 0.95
- Diferencias entre medias - Tabla:
tidy(TukeyHSD(anova_cp2))
- Diferencias entre medias - Gráfico:
tidy(TukeyHSD(anova_cp2)) %>%
ggplot(aes(x = contrast, y = estimate, ymin = conf.low, ymax = conf.high)) +
geom_point() +
geom_errorbar(width = 0.2) +
geom_hline(yintercept = 0, linetype = 2, color = "red") +
coord_flip()

par(mfrow = c(2, 2))
plot(anova_cp2)

Anova CP3
- Resultado análisis de varianza:
anova_cp3 <- aov(cp3 ~ Moments, data = datos)
tidy(anova_cp3)
emmeans(anova_cp3, specs = "Moments")
## Moments emmean SE df lower.CL upper.CL
## Habituation -0.2515 0.240 281 -0.7231 0.2201
## Treatment 0.1786 0.103 281 -0.0243 0.3816
## Break -0.1892 0.120 281 -0.4250 0.0466
## Final 0.0438 0.240 281 -0.4278 0.5154
##
## Confidence level used: 0.95
- Diferencias entre medias - Tabla:
tidy(TukeyHSD(anova_cp3))
- Diferencias entre medias - Gráfico:
tidy(TukeyHSD(anova_cp3)) %>%
ggplot(aes(x = contrast, y = estimate, ymin = conf.low, ymax = conf.high)) +
geom_point() +
geom_errorbar(width = 0.2) +
geom_hline(yintercept = 0, linetype = 2, color = "red") +
coord_flip()

par(mfrow = c(2, 2))
plot(anova_cp3)

Modelos Mixtos
Variación individual (1)
datos %>%
mutate(Individuo = as_factor(Individuo)) %>%
group_by(Moments, Individuo) %>%
summarise(promedio_cp = mean(cp1)) %>%
ggplot(aes(x = Moments, y = promedio_cp, color = Individuo, group = Individuo)) +
geom_point() +
geom_line() +
labs(color = "") +
guides(colour = guide_legend(nrow = 3)) +
theme(legend.position = "top")

datos %>%
mutate(Individuo = as_factor(Individuo)) %>%
group_by(Moments, Individuo) %>%
summarise(promedio_cp = mean(cp2)) %>%
ggplot(aes(x = Moments, y = promedio_cp, color = Individuo, group = Individuo)) +
geom_point() +
geom_line() +
labs(color = "") +
guides(colour = guide_legend(nrow = 3)) +
theme(legend.position = "top")

Variación individual (2)
datos %>%
mutate(Individuo = as_factor(Individuo)) %>%
group_by(Moments, Individuo) %>%
summarise(promedio_cp = mean(cp1)) %>%
ggplot(aes(x = Moments, y = promedio_cp, color = Individuo, group = Individuo)) +
facet_wrap(~Individuo, scales = "free", ncol = 5) +
geom_point() +
geom_line() +
labs(color = "") +
theme(legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1))

datos %>%
mutate(Individuo = as_factor(Individuo)) %>%
group_by(Moments, Individuo) %>%
summarise(promedio_cp = mean(cp2)) %>%
ggplot(aes(x = Moments, y = promedio_cp, color = Individuo, group = Individuo)) +
facet_wrap(~Individuo, scales = "free", ncol = 5) +
geom_point() +
geom_line() +
labs(color = "") +
theme(legend.position = "none",
axis.text.x = element_text(angle = 45, hjust = 1))

datos %>% filter(Pieza == "Pieza4")
Modelo CP1
library(nlme)
mod_cp1 <- lme(cp1 ~ Moments, method = "REML",
random = ~ Moments | Individuo,
data = datos,
control = lmeControl(opt = "optim"))
tidy(car::Anova(mod_cp1, type = "II"))
emmeans(mod_cp1, specs = "Moments")
## Moments emmean SE df lower.CL upper.CL
## Habituation -1.06 0.453 29 -1.989 -0.135
## Treatment 1.22 0.313 29 0.582 1.864
## Break -0.97 0.232 29 -1.446 -0.495
## Final -1.99 0.453 29 -2.913 -1.059
##
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
- Diferencias entre medias - Tabla:
tidy(pairs(emmeans(mod_cp1, ~Moments))) %>%
dplyr::select(contrast, estimate, std.error, adj.p.value)
- Diferencias entre medias - Gráfico:
library(multcomp)
library(multcompView)
medias <- as.data.frame(cld(emmeans(mod_cp1, ~Moments),
alpha = 0.05, Letters = letters,
adjust = "tukey"))
medias %>%
ggplot(aes(x = Moments, y = emmean, ymin = emmean - SE,
ymax = emmean + SE)) +
geom_point(position = position_dodge(0.9), alpha = 1) +
geom_errorbar(width = 0.1) +
geom_line(aes(group = 1)) +
geom_text(aes(label = .group, y = emmean),
position = position_dodge(0.9),
show.legend = FALSE, color = "black",
hjust = -0.3,
vjust = -0.5) +
scale_color_grey() +
scale_fill_grey() +
geom_hline(yintercept = 0, color = "black", linetype = 2) +
labs(y = "", x = "Moments") +
theme_bw() +
theme(legend.position = "none")

par(mfrow = c(1, 2))
#plot(mod_cp1)
residuales <- mod_cp1$residuals[, 1]
ajustados <- mod_cp1$fitted[, 1]
qqnorm(residuales)
qqline(residuales)
plot(x = ajustados, y = residuales)
abline(h = 0, col = "red", lty = 2)

Modelo CP2
mod_cp2 <- lme(cp2 ~ Moments, method = "REML",
random = ~ Moments | Individuo,
data = datos,
control = lmeControl(opt = "optim"))
tidy(car::Anova(mod_cp2, type = "II"))
emmeans(mod_cp2, specs = "Moments")
## Moments emmean SE df lower.CL upper.CL
## Habituation -1.009 0.316 29 -1.657 -0.362
## Treatment 0.967 0.185 29 0.588 1.347
## Break -0.657 0.158 29 -0.981 -0.334
## Final -1.288 0.323 29 -1.948 -0.627
##
## Degrees-of-freedom method: containment
## Confidence level used: 0.95
- Diferencias entre medias - Tabla:
tidy(pairs(emmeans(mod_cp2, ~Moments))) %>%
dplyr::select(contrast, estimate, std.error, adj.p.value)
- Diferencias entre medias - Gráfico:
medias2 <- as.data.frame(cld(emmeans(mod_cp2, ~Moments),
alpha = 0.05, Letters = letters,
adjust = "tukey"))
medias2 %>%
ggplot(aes(x = Moments, y = emmean, ymin = emmean - SE,
ymax = emmean + SE)) +
geom_point(position = position_dodge(0.9), alpha = 1) +
geom_errorbar(width = 0.1) +
geom_line(aes(group = 1)) +
geom_text(aes(label = .group, y = emmean),
position = position_dodge(0.9),
show.legend = FALSE, color = "black",
hjust = -0.3,
vjust = -0.5) +
scale_color_grey() +
scale_fill_grey() +
geom_hline(yintercept = 0, color = "black", linetype = 2) +
labs(y = "", x = "Moments") +
theme_bw() +
theme(legend.position = "none")

par(mfrow = c(1, 2))
#plot(mod_cp1)
residuales2 <- mod_cp2$residuals[, 1]
ajustados2 <- mod_cp2$fitted[, 1]
qqnorm(residuales2)
qqline(residuales2)
plot(x = ajustados2, y = residuales2)
abline(h = 0, col = "red", lty = 2)

Exportando gráficos
# Exportando gráfico pdf
ggsave(filename = "graphics_paper/acp.pdf",
plot = g1,
device = "pdf",
dpi = 96,
units = "cm",
width = 35.14, # igual a 531px
height = 14.05) # igual a 1328px
# Exportando gráfico Tiff
ggsave(filename = "graphics_paper/acp.tiff",
plot = g1,
device = "tiff",
dpi = 96,
units = "cm",
width = 35.14, # igual a 531px
height = 14.05) # igual a 1328px